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Watershed  Scale  TMDL  Model:  Multidimensional 
Sediment  Erosion,  Transport,  and  Fate 


By  Charles  W.  Downer  and  Aaron  Byrd 


PURPOSE:  The  purpose  of  this  System-Wide  Water  Resources  Program  (SWWRP)  technical 
note  is  to  describe  the  new  sediment  transport  routines  in  GSSHA,  a  watershed  analysis  and 
management  tool  that  has  the  ability  to  simulate  the  movement  of  water,  sediment,  and 
associated  constituents  at  fine-scale  increments  (<100  m)  over  watershed  scale  areas.  The 
resulting  tool  is  intended  for  analyzing  project  alternatives  and  best  management  practices 
(BMPs)  designed  to  control  sediments  near  the  source,  either  on  upland  areas  or  in  tributaries. 

BACKGROUND:  Receiving  water  bodies  are  harmed  by  the  introduction  of  excess  sediments 
which  reduce  storage  capacity,  increase  turbidity,  and  introduce  associated  contaminants.  The 
control  of  sediments  may  best  be  performed  in  upland  areas  near  their  source  where 
contaminants  may  be  removed  along  with  the  sediments  before  the  contaminants  are  released 
into  solution. 

In  this  work  unit,  the  Department  of  Defense  Watershed  Modeling  System  (WMS)-distributed 
hydrologic  model  Gridded  Surface  Subsurface  Hydrologic  Model  (GSSHA)  (Downer  et  al. 

2005)  has  been  further  developed  to  allow  the  physics  based  simulation  of  sediment  erosion, 
transport,  and  deposition  on  a  continuous  basis.  Simulation  of  sediments  is  at  the  most 
fundamental  level  of  current  understanding  of  sediment  physics.  The  general  laws  of 
conservation  of  mass  and  momentum  are  applied.  This  approach  allows  contaminants  closely 
associated  with  sediments  to  be  simulated  with  the  same  approach  and  equations.  A  continuous 
simulation  model  provides  the  ability  to  track  sediments  over  several  precipitation  events,  or 
years,  and  determine  the  long-tenn  patterns  of  erosion,  deposition,  and  geomorphologic  changes 
within  a  watershed. 

GSSHA  is  a  distributed,  physics-based  hydrologic  model  developed  to  simulate  a  watershed’s 
response  to  meteorological  inputs.  Basic  simulated  physical  processes  include  distributed 
rainfall,  rainfall  interception  by  vegetation,  surface  ponding  and  retention,  infiltration, 
evapotranspiration,  overland  flow,  streamflow,  and  lateral  saturated  groundwater  flow.  The 
model  also  simulates  subsurface  drainage  networks  and  includes  seasonality  effects.  The  model 
is  intended  for,  and  has  been  applied  to,  the  simulation  of  streamflow,  flooding,  soil  moisture, 
sediment  erosion,  and  discharge. 

The  hydrologic  and  hydraulic  components  of  GSSHA  provide  the  information  necessary  to 
simulate  erosion  and  sediment  fate  and  transport.  Key  to  simulating  sediments  within  a 
watershed,  GSSHA  provides  the  ability  to  link  two-dimensional  (2-D)  overland  flow  transport  to 
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one-dimensional  (1-D)  transport  in  a  stream  network  as  shown  in  Figure  1.  The  overland  flow 
plane  provides  inputs  to  the  channel  network  at  each  point  where  the  channel  network  and 
overland  flow  plane  coincide.  Optionally,  water  and  constituents  can  be  allowed  to  spill  back  on 
the  overland  flow  plane. 


Flow  in  both  domains  is  described  by  the  diffusive  wave  approximation  of  the  de  Saint  Venant 
equations  and  solved  using  the  finite  volume  technique.  The  finite  volume  method  allows  the 
simultaneous  simulation  of  both  wet  and  dry  cells.  The  ability  to  simulate  “dry  bed”  conditions 
is  critical  in  the  simulation  of  large  regional  watersheds  as  precipitation  may  occur  preferentially 
within  the  study  area  and  only  parts  of  the  watershed,  or  stream  network,  may  be  flowing  at  any 
given  time. 

Prior  to  this  effort,  GSSFIA  employed  an  empirical  approach  to  sediment  erosion,  transport,  and 
fate  (Johnson  et  al.  2000;  Sanchez  2002)  taken  from  the  CASC2D  model  (Ogden  and  Julien 
2002),  GSSFIA’s  predecessor.  Due  to  limitations  of  this  method  (Ogden  and  Heilig  2001)  and  a 
desire  to  simulate  sediments  and  dissolved  constituents  with  one  general  method,  the  more 
empirical  method  has  been  replaced  with  a  more  general,  physics-based  approach  as  described  in 
this  document. 


2 


ERDC  TN-SWWRP-07-1 
January  2007 


METHODOLOGY:  The  general  transport  equations  are  solved  in  two  dimensions  on  the 
overland  flow  plane  and  in  one  dimension  in  the  channel  network,  using  a  common  methodology 
and  interaction  between  the  two  domains.  Sediment  erosion  and  deposition  are  treated  as  first 
order  reactions,  with  erosion  being  an  source  tenn,  and  deposition  being  a  sink  term. 

GENERAL  TRANSPORT  EQUATIONS:  The  general  transport  equations  describe  the  fate 
and  transport  of  dissolved  and  particulate  constituents.  Each  cell  is  treated  as  a  completely 
mixed  reactor  and  the  concentration  in  each  cell  is  affected  by  internal  sources  and  sinks,  as  well 
by  advection  into  and  out  of  the  cell  and  diffusion  between  surrounding  cells. 

For  the  sake  of  simplicity,  the  1-D  scheme  (channels)  is  presented  first.  Although  only  one 
constituent  is  described  in  this  text,  the  number  of  constituents,  both  dissolved  and  particulate 
(sediment),  can  be  any  number  specified  by  the  user. 

Reactive  Transport  in  Channels.  GSSHA  computes  values  of  water  flow  and  depth  within 
a  user-specified  1-D  finite  volume  stream  network.  Solution  of  the  diffusive  wave 
approximation  provide  the  discharge  and  cross-sectional  area  at  fine  space  and  time  increments 
within  the  channel  network.  Typical  stream  nodes  range  in  size  from  30  to  200  m.  Typical 
channel  routing  time-steps  range  from  5  to  30  sec. 


This  information  is  needed  for  the  transport  of  sediments  and  other  constituents  within  the 
channel  network  using  the  general  1-D  advection-dispersion  equation  in  terms  of  the  mass  of 
constituent  (M)  equal  to  the  concentration  (C)  multiplied  by  the  volume  ( V)  with  constant 
dispersion: 


d(M)  |  d(uM) 
dt  dx 


—fn— 

dx  v  dx  j 


+  K(M)  =  S 


(1) 


where: 


u  =  flow  velocity  in  the  x-direction  (M  T'1) 

C  =  concentration  (M  L'3) 

Qx  =  discharge  in  the  x-direction  (L  T’  ) 

A  =  cross-sectional  area  (L“) 

Dx  =  diffusion  coefficient  in  the  x-direction  for  the  constituent  of  concern  (L  T"  ) 
K  =  decay  coefficient  (T1) 

V  =  volume  (L3) 

S  =  term  for  all  sources  and  sinks  (M  T'1) 


This  equation  is  solved  with  a  simple  three-point  explicit  finite -volume  scheme  for  each 
constituent  of  concern.  The  solution  proceeds  from  the  most  upstream  node  and  proceeds 
downstream  to  the  channel  outlet.  This  scheme  is  first  order  accurate  in  time  and  space.  For 
each  node,  the  following  mass  reactions  are  accounted  for  with  masses  ( M)  in  grams;  volumes 
(V)  in  m  ;  discharges  ( Q )  in  m  s’  ;  and  concentrations  (C)  in  g  m"  or  mg  L‘  . 
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•  Initial  mass, Ma  =  Ax(C"A") 

o  Where  A  is  the  cross-sectional  area,  Ax  is  the  channel  grid  element  size,  the 

superscript  n  refers  to  the  current  time  level,  and  the  subscript  i  refers  to  the  current 
node. 

•  Advective  flux,  MaJ,„  =  Affec”,  -  Q"'C 

o  Where  At  is  the  time-step  (s);  the  subscripts  i -1,  and  i  refer  to  upstream  node,  and 
current  node,  respectively,  with  flow  being  defined  as  positive  from  the  current  cell  to 
the  downstream  cell;  the  superscript  n+/  refers  to  the  next  time  level. 

•  Constant  point  source,  M p  =  A t(QpC p ) 

o  Where  Qp  and  Cp  are  the  flow  and  concentration  of  the  point  source,  respectively. 

•  Lateral  flow,  exchange  between  channel  and  overland  flow  plane,  Mlat  =  AtAl(qJ+lC"at ) . 

o  Where  A /  is  length  of  the  stream/overland  interface  (m),  and  Ciat  is  the  concentration 
of  the  water  on  the  overland  plane  if  flow  is  into  the  channel,  and  the  concentration  of 
the  water  in  the  channel  if  flow  is  back  onto  the  overland  plane. 

•  Groundwater  exchange,  M ^  =  At(g"+'C”1V) 

o  Where  the  concentration  of  the  groundwater,  Cgw,  depends  on  the  direction  of  Q.  It  is 
equal  to  C,  if  the  flux  is  into  the  groundwater  (Q  is  negative)  or  equal  to  concentration 
in  the  groundwater  if  the  flux  is  into  the  stream  (Q  is  positive). 


Dispersive  exchange, =^(A.„r4"',2(CA, -C')- D^A^C"  -C’,)) 


Ax 


9  1 

o  Where:  D  is  the  dispersion  coefficient  (m  s'  ),  assuming  D  varies  in  space,  but  not  in 
time,  and  Ai+i/2  and  Aui/2  are  the  cross-sectional  areas  (m2)  between  the  i  and  i+1  and 
i  and  i-1  nodes,  respectively. 


•  Decay,  Mfc,=A/(A-"C;T") 

o  Where  Kdecay  is  the  decay  coefficient  computed  at  each  time  increment  based  on  the 
physical  state  of  the  sytem  (s'1). 


Concentration  at  the  next  time  level,  n+1,  is  calculated  as: 


c;+1 


Mo  +Madvecl  +  Mp  +Mlat  +  M  gw  +  M  disp  +M  decay 

V"H 


(2) 


The  mass  exchange  between  the  channel  and  the  groundwater  is  only  computed  if  saturated 
groundwater  computations  are  being  conducted.  The  concentration  in  the  groundwater  is 
assumed  to  be  zero,  so  that  Cgw  =  0  when  the  stream  is  gaining.  The  concentration  of  the 
overland  flow  plane,  Ciat,  will  also  be  zero  unless  the  overland  flow  transport  is  also  being 
simulated. 
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Overland  Reactive  Transport.  For  overland  flow,  the  2-D  form  of  the  general  transport 
equation  is  solved. 


d(M)  d(uM)  d(yM)  _  d  f  dM ' 
dt  dx  dy  dx  v,  r  dx  , 


d_ 

dy 


D, 


dM ' 
dy  y 


+  KM  = 


5 


(3) 


Where  v  is  the  velocity  in  the  y-direction  (m  s’1). 


A  2-D  overland  flow  solution  of  the  diffusive  wave  approximations,  as  described  by  Downer  and 
Ogden  (2004),  provides  the  necessary  discharges  and  areas.  As  with  the  instream  contaminant 
transport  routines,  the  overland  flow  transport  is  modeled  with  a  mass  balance  in  each  overland 
flow  cell.  For  each  constituent  and  each  cell,  the  following  mass  reactions  are  accounted  for 
with  masses  (M)  in  grams;  volumes  (  V)  in  nr  ;  discharges  (Q)  in  m  s’  ;  and  concentrations  (C) 
in  g  m’3  or  mg  L’1. 


•  Initial  mass, Mn  =  AtAx1C"h"i 

o  J  U 

o  Where:  Ax  is  the  grid  size  (m),  h  is  depth,  and  the  subscript  ij  refers  x  and  y  location, 
respectively,  of  the  current  node. 

•  Precipitation  Mass,  M pre  =  AtAx2  (IC pre) 

o  Where  /  is  rainfall  intensity  (m  s’1)  and  Cpre  is  the  concentration  in  the  precipitation 
(assumed  to  be  zero). 

•  Advective  flux  in  the  x-direction,  Mx  =  AtAx(p.^'JC',_lj  -  p-jC-j ) 
o  Where  p  is  the  unit  discharge  (uh)  in  the  x-direction  (m  s’  ). 

•  Advective  flux  in  the  y-direction,  M y  =  A/Ax(y''JllC("/  ,  -  q"+.  C;" ) 
o  Where:  q  is  the  unit  flux  (yh)  in  the  y-direction  (m  s’  ). 

•  Point  source,  Mp  -  At (<2pCp), 


•  Lateral  inflow  to  channel,  Mlat  =  -AtAl[q"^C"at) 

•  Infdtration  loss,  Minf  -  AtAx2  (/„"+1C” ) 

o  Where:  f„f  is  the  infiltration  rate  (m  s’1). 

•  Exfiltration,  M exf  =  AtAx2(f^'Cgw) 

o  Where:  fexf  is  the  exfiltration  rate  (m  s’1)  and  Cgw  is  assumed  to  be  constant. 


Dispersive  exchange  in  x-direction 
At 


M 


dispx 


Ax' 


\d 


jn+\ 

f+l/2Jyif+l/2J 
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•  Dispersive  exchange  in  the  y-direction, 

-cp-D^A^c;  -c;y,)) 

•  Decay,  Mimy=At(M,K;^) 


•  Uptake  from  the  land  surface,  M  uptake  =  K^ake  (Cmax  -  C” )  At  Ax 2 

2  1 

o  Where:  K uptake  is  the  uptake  coefficient  (g  m'  s'  )  and  Cmax  is  the  maximum  possible 
concentration  in  the  overland  flow  cell. 

The  2-D  overland  transport  equations  are  solved  with  a  five-point  explicit  alternating  direction 
(ADE)  scheme  with  prediction  correction  (PC),  or  ADE-PC  (MacCormick  1971).  This  method 
is  2nd  order  accurate  in  space  and  time.  The  concentration  at  the  next  time-step,  n+1  time  level,  is 
computed  using  the  following  method. 


•  If  the  overland  cell  is  located  along  a  channel,  the  volume,  mass,  and  concentration  in 
that  cell  is  first  corrected  for  lateral  inflow  into  the  channel  node  or  channel  flow  back  on 
the  overland  flow  plane. 

•  The  volume,  mass,  and  concentration  are  then  adjusted  for  precipitation,  point  sources, 
and  exfiltration. 

•  Uptake  and  decay  rates  are  calculated  in  each  node,  solids  only,  as  described  in  the  next 
section. 

•  Infiltration,  dispersion,  decay  and  advection  in  the  x-  and  y-directions  are  computed  using 
an  alternating  direction  explicit  ADE-PC,  which  proceeds  in  the  following  manner. 

o  Changes  in  cell  volume,  mass  and  concentration  due  to  infiltration,  decay,  and 
advection  and  dispersion  in  the  y-direction  are  first  computed  using  values  at  the  n 
time  level. 

o  Intermediate  values  of  volume  and  concentration  at  the  n+1/2  time  level  are 
computed  based  on  the  mass  changes  at  the  n  time  level, 
o  An  average  of  the  values  of  volumes  and  concentrations  at  the  n  time  level  and  the 
intermediate  values  at  the  n+1/2  time  level  are  used  to  compute  new  values  of 
concentration  at  the  n+1/2  time  level. 

o  Changes  in  mass  and  concentration  due  to  infiltration,  decay,  and  advection  and 
dispersion  in  the  y-direction  are  computed  using  the  new  values  of  concentration  at 
the  n+1/2  time  level. 

o  Final  values  of  concentration  at  the  n+1/2  time  level  are  computed  based  on  the  final 
mass  fluxes  computed  at  the  n+1/2  time  level, 
o  The  preceding  steps  are  repeated  for  the  x-direction,  yielding  values  of  concentration 
at  the  n+1  time  level.  To  avoid  any  directional  bias,  the  order  of  solution,  first  v  and 
thenx,  or  first  x  and  then y,  changes  every  time-step. 
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When  solids  are  being  transported,  there  is  no  loss  or  gain  of  the  solids  as  water  infiltrates, 
evaporates,  rains,  or  exfiltrates,  but  the  concentration  of  solids  in  the  cell  is  modified  in  response 
to  the  reduction  in  the  volume  of  water  in  the  cell. 

Uptake  —  Erosion.  Though  the  physics  of  erosion  on  the  overland  flow  plane  are  not  well 
described,  noncohesive  particles  are  thought  to  erode  due  to  the  hydraulic  properties  of  flow  and 
rainfall  impact.  However,  in  order  to  test  the  overall  model  formulation,  some  currently  applied 
methods  of  computing  overland  erosion  are  incorporated  into  the  model. 

In  overland  flow  simulations,  the  erosion  is  typically  described  in  tenns  of  the  transport  capacity, 
the  amount  of  sediment  that  can  be  transported  under  the  given  hydraulic  conditions.  As 
reviewed  by  Julien  and  Simons  (1985),  transport  capacities  are  typically  of  the  form  of  the 
original  Kilinc -Richardson  equation  (Kilinc  and  Richardson  1973): 

qs  =  25500g2035^'664  (4) 

where  the  factor  25,500  is  an  empirical  constant,  q,  is  the  sediment  unit  discharge  (ton  m'1  s' 1 )  q 
2  1 

is  unit  discharge  (nr  s'  )  and  S0  is  the  land  surface  slope. 

Julien  (1995)  modified  the  original  Kilinc-Richardson  equation  to  allow  the  transport  capacity  to 
be  computed  with  conditions  of  nonunifonn  flow  with  consideration  of  soil  and  land-use  specific 
factors  to  compute  the  transport  capacity  (tons  m"1  s'1): 

qH  =25500 q2m5S\6M  E*C*P  (5) 

f  0.15 

where: 

q  =  unit  discharge  (m  s'  ) 

Sf  =  friction  slope  (dimensionless) 

E  =  soil  erodability  factor,  with  values  ranging  from  0  to  1 
C  =  soil  cropping  factor  (0-1) 

P  =  conservation  factor  (0-1) 

In  CASC2D-SED  (Johnson  et  al.  2000;  Sanchez  2002)  and  previous  versions  of  GSSHA,  this 
transport  capacity  is  computed  in  both  the  x-  and  y-directions  and  then  used  to  transport 
deposited  and  parent  materials,  proportional  to  available,  in  each  direction.  The  transport 
capacity  is  always  satisfied. 

For  inclusion  into  the  preceding  general  transport  equations,  the  transport  capacity  must  be 
changed  into  an  erosion  coefficient.  Because  determining  the  interdependent  factors  E,  C,  and  P 
is  difficult  (Ogden  and  Helig  2001),  they  are  combined  into  a  single  coefficient,  Fero<je, 
representing  the  overall  erodability  of  soils.  The  x  and  y  components  of  both  unit  discharge,  q, 
and  friction  slope,  Sj,  are  used  to  compute  the  magnitude  of  total  transport  capacity  in  the  Kilinc- 
Richardson  equation.  When  divided  by  the  grid  size,  Ax,  an  uptake  rate,  Kuptake,  (g  m'2  s'1)  is 
provided  to  be  used  in  the  general  transport  equations.  The  final  equation  is: 
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K 


uptake 


25500 
1. 102x1  (T6 


2.035  c  1.664 

q  Sf 


F 

erode 

0.15Ax 


(6) 


Alternately  the  uptake  rate  for  each  particle  class  being  simulated  may  be  computed  from  the 
Engelund-Hansen  equation  (Engelund  and  Hansen  1967)  as: 


K 


{)^Ferodeysv 


uptake 


Ax 


d 


50 


g(^-D 

r, 


Ml 


-i3/2 


(7s  ~  Yf  Ko 


where 

v 

dso 

g 

Yf 

Y 

To 


magnitude  of  flow  velocity  (m  s'1) 

median  particle  size  for  each  class  specified  (m) 

acceleration  of  gravity  (m  s'2) 

-5 

specific  weight  of  the  fluid  (N  m"  ) 
specific  weight  of  the  particle  (N  m"  ) 
bed  shear  stress  (N  m'  ),  computed  as: 


T0  =  YfhSf 


(7) 


(8) 


where:  h  is  the  depth  of  water  in  the  cell  of  interest  (m)  which  may  range  from  less  than  a 
centimeter  in  areas  of  sheet  flow  to  perhaps  one  meter  in  concentrated  flow  areas.  The  velocity 
in  the  cell  can  be  computed  from  the  x  and  y  components  of  unit  discharge  as: 


\(  \ 

2 

r  „  \ 

[  p 

+ 

Jl  j 

1  h  J 

(9) 


According  to  Govers  (1990),  transport  does  not  occur  unless  the  stream  power,  Q,  is  greater  than 
0.004  m  s'1.  Therefore,  for  the  Engelund-Hansen  equation,  Kuptai \e  =  0  for  Q  <  0.004.  The  stream 
power  is  computed  as: 

Q  =  vSf  (10) 

Erosion  due  to  rainfall  impact  is  simulated  as  described  by  Dario  (2002).  When  total  erosion  is 
calculated,  as  in  the  Kilinc-Richardson  equation,  the  amount  of  erosion  of  each  particle  size  is 
determined  from  the  fractions  of  particles  available  in  the  parent  material.  For  cells  with 
deposited  materials,  these  materials  are  eroded  first,  then  any  additional  erosion  that  occurs  is 
from  the  parent  material.  The  distribution  of  particles  in  the  deposited  materials  is  likely 
different,  and  tracked  separately,  from  the  parent  material.  Particles  are  tracked  according  to  size 
(as  specified),  location  (which  cell),  and  status  (parent  material,  in  suspension,  or  deposited). 
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Maximum  Concentration.  Uptake  mass  is  limited  by  the  maximum  concentration  that  the 
fluid  can  transport.  This  can  be  taken  as  a  concentration  equal  to  the  specific  gravity  of  the 
particle  times  the  water  density,  or  computed  from  a  variety  of  methods.  As  implemented  in 
KINEROS  (Woolhiser,  Smith,  and  Goodrich  1990),  the  Engelund-Hansen  equation  can  be  used 
to  determine  the  maximum  concentration,  Cmax,  (g  in"  )  for  each  sediment  as: 


C 


max 


0.05SGp 
d5o  (SG  - 1)2 


0.004) 


(ID 


where  SG  is  the  specific  gravity  of  the  particles,  and  p  is  the  density  of  water  at  20°  C  (g  m"3). 

These  equations  are  included  because  they  have  some  general  applicability  and  acceptance  and 
some  method  of  estimating  erosion  uptake  rate  and  maximum  concentration  is  needed  to  allow 
the  overall  method  to  be  tested. 


Decay  -  Deposition.  Deposition  is  controlled  by  particle  shape,  size  (dso),  specific  gravity 
(SG),  and  the  properties  of  the  fluid  that  the  particle  is  settling  in.  For  particles  larger  than  1  mm, 
the  settling  velocity,  vs  (m  s"1)  is  (Julien  1995): 


co- 


8v 

d 


yjl  +  0.0139J*  -1 


(12) 


where 

2  1 

v  =  kinematic  viscosity  of  the  fluid  (m  s"  ) 
d  =  median  particle  size,  d50  (m) 
d*  =  dimensionless  particle  diameter,  computed  as: 


(SG-\)g 

,,2 


(13) 


For  particles  smaller  than  1  mm,  the  Stokes  equation  is  used: 


co  = 


gd2  7s -if 

i8v  rf 


(14) 


The  kinematic  viscosity  of  water  depends  on  water  temperature.  For  long-tenn  simulations, 
where  hourly  values  of  air  temperature  are  provided,  overland  water  temperature  is  assumed  to 
be  equal  to  the  air  temperature.  For  single  event  simulations  the  water  temperature  is  specified; 
20°  C  is  the  default  value. 


In  the  general  transport  equations  previously  described,  deposition  is  treated  as  decay.  The 
decay  rate  Kdecay  (s  ')  is  calculated  each  time-step  as: 
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Kdecay=^  (15) 

REACTIVE  TRANSPORT:  Although  the  schemes  described  in  this  technical  note  were 
developed  for  transport  of  sediments,  the  techniques  employed  are  general  and  the  same  methods 
could  be  used  to  simulate  the  fate  and  transport  of  dissolved  reactive  constituents.  For  simple 
first  order  reactions,  the  decay  rates  ( K)  would  reflect  decay  rates  for  the  loss  of  the  reactive 
constituent.  For  complex  higher  order  reactions  with  multiple  subspecies  and  transformations 
between  these  species,  the  change  in  mass  due  to  multiple  reactions  (KV)  could  be  calculated 
outside  the  GSSHA  model  and  passed  back  to  the  reactive  transport  routine.  This  approach 
would  allow  any  number  of  subspecies  to  be  accounted  for  without  adding  undue  complexity  to 
the  basic  GSSHA  formulation. 

SUMMARY :  The  GSSHA  model  has  been  modified  to  allow  the  transport  and  fate  of 
constituents  on  both  the  overland  flow  plane  and  within  the  stream  network  with  the  general 
transport  equations.  The  user  specifies  the  number  and  properties  of  constituents  to  be 
simulated.  The  uptake  of  constituents  from  the  overland  flow  plane  and  loss  of  materials  in 
either  the  overland  flow  plane  or  stream  network  is  simulated  as  first  order  reactions.  For  solids, 
sediments,  erosion  by  any  means,  hydraulic  or  due  to  rainfall  impact,  is  converted  into  an  uptake 
rate.  Deposition  is  calculated  as  a  decay  rate.  These  uptake  and  decay  rates  are  computed  in 
each  grid,  during  each  time-step.  The  method  allows  sediments  to  erode  in  one  cell,  be  advected 
to  downstream  cells,  and  be  deposited  in  either  the  cell  of  origin  or  in  downstream  cells.  Erosion 
can  be  simulated  for  extended  periods,  allowing  the  trends  in  erosion  and  deposition,  and 
resulting  morphological  changes  in  the  watershed  to  be  tracked.  Ongoing  research  will  identify 
the  important  processes  in  erosion,  transport,  and  deposition  providing  improved  estimates  of 
erosion  rates  and  sediment  redistribution  on  the  landscape. 

ADDITIONAL  INFORMATION:  This  SWWRP  technical  note  was  prepared  by  Dr.  Charles  W. 
Downer  and  Aaron  Byrd,  Coastal  and  Hydraulics  Laboratory,  U.S.  Army  Engineer  Research  and 
Development  Center.  The  study  was  conducted  as  an  activity  of  the  Regional  Sediment 
Management  work  unit  of  the  System-Wide  Water  Resources  Program  (SWWRP).  For 
information  on  SWWRP,  please  consult  https://swwrp.usace.armv.mil/  or  contact  the  Program 
Manager,  Dr.  Steven  L.  Ashby  at  Steven.L.Ashbv@erdc.usace.army.mil.  This  technical  note 
should  be  cited  as  follows: 

Downer,  C.  W.,  and  A.  R.  Byrd.  2007.  Watershed  scale  TMDL  model:  Multi¬ 
dimensional  sediment  erosion,  transport,  and  fate.  ERDC  TN-SWWRP-07-1. 

Vicksburg,  MS:  U.S.  Army  Engineer  Research  and  Development  Center. 

https://swwrp.usace.army.mil/ 
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